﻿#ifndef LDUDILUPRECOND_HPP
#define LDUDILUPRECOND_HPP

#include "unapMatrix.hpp"
#include "unapPreconditioner.hpp"

namespace UNAP {
template <typename Cmpt>
class LduMatrix;
template <typename Cmpt>
class CSRMatrix;

/// \brief Simplified Diagonal-based Incomplete LU preconditioner
/// support symmetric and asymmetric, same as DIC if matrix is symmetric
template <typename Cmpt>
class DILUPrecond : public Preconditioner<Cmpt> {
 private:
  /// \brief the reciprocal diagonal
  Vector<Cmpt> rD_;

  /// \brief LduMatrix info
  label nFaces_;
  label *uPtr_;      /*!< upper addr */
  label *lPtr_;      /*!< lower addr */
  label *losortPtr_; /*!< lower sort addr */
  Cmpt *upperPtr_;   /*!< upper value */
  Cmpt *lowerPtr_;   /*!< lower value */
  bool ownData_;

 public:
  /// \brief constructor
  /// \param A matrix input
  DILUPrecond(const LduMatrix<Cmpt> &A);
  DILUPrecond(const CSRMatrix<Cmpt> &A);

  /// \brief destructor
  virtual ~DILUPrecond() {
    if (ownData_) {
      DELETE_POINTERS(uPtr_)
      DELETE_POINTERS(lPtr_)
      DELETE_POINTERS(upperPtr_)
      DELETE_POINTERS(lowerPtr_)
      DELETE_POINTERS(losortPtr_)
    }
  }

  /// \brief calculate the reciprocal of the preconditioned diagonal
  /// \param rD the resulting reciprocal diagonal
  /// \param A matrix input
  void calcReciprocalD(Vector<Cmpt> &rD);

  /// \brief return wA the preconditioned form of residual rA
  /// \param w resulting vector, preconditioned by rA
  /// \param r Preconditioner vector
  virtual void precondition(Vector<Cmpt> &w, const Vector<Cmpt> &r) const;

  /// \brief return the reciprocal diagonal
  virtual const Vector<Cmpt> &rD() const { return rD_; }
};

}  // namespace UNAP

#endif  // LDUDILUPRECOND_HPP
